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! ! Abstract 

(— I ', 

A quantum control landscape is defined as the expectation value of a target 
j \ | observable as a function of the control variables. In this work control landscapes 

for open quantum systems governed by Kraus map evolution are analyzed. Kraus 
maps are used as the controls transforming an initial density matrix p x into a final 
density matrix to maximize the expectation value of the observable 0. The absence 
of suboptimal local maxima for the relevant control landscapes is numerically illus- 
I ■ trated. The dependence of the optimization search effort is analyzed in terms of the 

dimension of the system N, the initial state p\, and the target observable 0. It is 
■^■j- \ found that if the number of nonzero eigenvalues in p x remains constant, the search 

effort does not exhibit any significant dependence on N. If p\ has no zero eigenval- 
ues, then the computational complexity and the required search effort rise with N. 
The dimension of the top manifold (i.e., the set of Kraus operators that maximizes 
' the objective) is found to positively correlate with the optimization search efficiency 

<^ • Under the assumption of full controllability, incoherent control modelled by Kraus 

J> \ maps is found to be more efficient in reaching the same value of the objective than 

V^j ' coherent control modelled by unitary maps. Numerical simulations are also per- 

formed for control landscapes with linear constraints on the available Kraus maps, 
and suboptimal maxima are not revealed for these landscapes. 



1 Introduction 

The general goal of quantum control is to apply a suitable external field to a system in 
order to maximize the expectation value of a target operator. If the system under control is 
isolated from the environment, then the dynamics are coherent and described by a unitary 
transformation, as appears in coherent control [U [2j [3j SI El [6], [TJ, [HI [9] - In practice, all 
real systems are open and interact with the environment in some fashion, so that the 
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dynamics of the system will have some incoherent component. Such incoherent control 
through interaction with the environment [TQ1 EH E21 ESI EU ESI ESI or through quantum 
measurements [171 ESI ESI EDI EH E21 EH] can be beneficial in many cases (e.g., in creation 
of thermal beams of metastable noble gases |24j . in quantum computing with mixed 
states [25], or in the modification [26] of Grover's algorithm to extend the capabilities of 
the original unitary scheme). 

In this paper, we consider the most general class of physically-allowed state trans- 
formations of controlled open quantum systems. These transformations are represented 
by Kraus maps [27J providing a kinematic description of incoherent control. Embedded 
in these maps is information about the system and environment, both of which may be 
subject to control. A control action determines the system's evolution with a Kraus map 
$, which transforms an initial system state p- x into the evolved final state pf = $(pi). The 
final state pf determines the expectation value J[$] := (6) = Tr[<3>(pi)0] = Tr[p f 0] of 
a target Hermitian operator representing a desired physical property to be optimized. 
The corresponding control goal is formulated as follows: given an initial state p x and a 
target observable 0, find a Kraus map $ opt that transforms p x into a state maximizing 

the expectation value, i.e., such that J[$ pt] = max «/[$]. The set of all Kraus maps for a 

■I' 

given quantum system forms a complex Stiefel manifold to formulate the control goal as 
a nonlinear problem of maximizing the objective function J over the Stiefel manifold. 

As shown in |28j . for any desired final state pf there exists a Kraus map <3> Pf that 
transforms all initial states pi into pf, i.e., such that $ Pf (p;) = Pf for any p y . If \ip) 
is an eigenvector of the target operator that corresponds to the maximal eigenvalue 
6U a x, then the expectation Tr[pf0] is maximized by the final state pf = pf^ = 
and therefore the objective J[$] = Tr[$(pi)0] is maximized, e.g., by the optimal map 
$ = $ Pf . The corresponding maximum objective value is J max = # max . Thus, the ability 
to generate dynamically arbitrary Kraus maps for an open quantum system implies its 
complete state-to-state controllability and, in particular, complete controllability for the 
objectives of the form J = (0). In contrast, a closed quantum system controlled by unitary 
dynamics has restricted state controllability; if pi and pf do not have the same eigenvalue 
spectrum, there does not exist a unitary transformation U such that JJ p\W = pf. The 
maximum attained value for the objective J = (0) in this case will generally be less than 

max- 

The quantum control landscape is defined by J = (0) as a function of the control vari- 
ables. The ability to successfully use a gradient or other local algorithm for maximization 
of the objective function depends on the existence or the absence of suboptimal local 
maxima. If local maxima exist, a local algorithm could get stuck at such points, and for 
this reason, we refer to suboptimal local maxima as "traps;"the presence of local saddle 
points should not serve as traps. In the case of coherent laser control, the landscape is 
known to be trap free [291 ES] . 

A detailed analysis of the control landscapes for incoherent control of open two-level 
quantum systems was performed [31] , where the absence of traps for these landscapes was 
proven. Arbitrary multi- level systems were considered in [32] , where it was shown that no 
suboptimal traps exist for the control landscapes for any finite-level open quantum system. 
In addition, a high-dimensional submanifold of optimal controls was found. As in the case 
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of coherent control, these results on the absence of traps and the multi-dimensionality of 
the global optimum manifold provide a theoretical foundation for the empirical fact that 
it is relatively easy to find optimal solutions even in the presence of an environment. 

The absence of traps in control landscapes for both closed and open quantum systems 
implies that the search using a local algorithm will eventually reach a global optimum 
solution. However, the absence of traps does not specify the efficiency of optimization 
procedure and the search effort needed to reach the solution. The efficiency of the op- 
timization procedure to find an optimal control, which is of practical importance due to 
limitations on computer time in simulations and laboratory resources in experiments, is 
determined by the local features of the control landscape as well as its topological char- 
acteristics. Different trap-free control landscapes can exhibit different degrees of search 
complexity. The prior relevant theoretical landscape analyses [3T| [32] for incoherent con- 
trol of open quantum systems did not describe the dependence of efficiency on the key 
parameters of the control problem: the dimension of the system N and the eigenval- 
ues of the target operator 6 and the initial state p. For closed quantum systems, a 
theoretical analysis of the computational complexity of coherent control landscapes was 
performed [331 US] along with a numerical analysis of the search effort using gradient, 
genetic and simplex algorithms [53, ES] • The results indicate that the search effort scales 
weakly, or possibly independently, with the dimension of the system N. 

This paper presents a numerical analysis, with a gradient algorithm, of the search 
effort for incoherent control of open quantum systems. The analysis lends insight into 
the topological and structural characteristics of the corresponding quantum control land- 
scapes. It shows that the search effort for driving a pure state into another pure state 
with Kraus maps remains relatively constant as the dimension N of the system increases, 
and this behaviour is qualitatively similar to the scaling behavior of the search effort for 
closed systems [35l [36] . A more general result is established for arbitrary, not necessarily 
pure, initial states: the search effort is essentially determined by the number of nonzero 
eigenvalues of the initial state p, and not by the dimension of the system N. Thus, when 
the number of non-zero eigenvalues of the initial state remains constant, the search effort 
does not depend on N. At the extreme of driving a mixed state with no zero eigenvalues 
into a pure state the search effort increases with the dimension of the system. The detailed 
analysis shows that the search effort is sensitive to the eigenstructure of the initial state p 
and the target operator 6; specifically, the degeneracies of the zero eigenvalue of p and of 
the maximal eigenvalue of positively correlate with the search efficiency, so that higher 
values of these degeneracies require less optimization search effort and correspond to a 
more efficient search. Further, comparative analysis of incoherent and coherent control 
shows that incoherent control under the full controllability assumption is a more efficient 
process than coherent control, indicating that the additional control freedom afforded by 
incoherent control can decrease the complexity of the problem. Finally, an analysis of 
control landscapes with linear constraints on the control variables is performed, and it 
does not reveal the presence of suboptimal traps even for a large number of independent 
constraints. Use of Kraus maps for modelling the controlled evolution of the system in 
this paper greatly simplifies computations as it does not require solving the dynamical 
evolution equations. Analysis of the scaling properties of the search effort for dynamical 
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optimization of open quantum systems remains as an issue for future study that can be 
performed using various specific models for the system and the environment [37J [55] . 

The paper is organized as follows. Section [2] describes the general theoretical frame- 
work for the kinematic analysis of incoherent control of multilevel open quantum systems. 
The expressions for the gradient and Hessian of the objective function J are derived in 
Sec. [3j Section H] contains the results of the numerical simulations. Section 4.1 describes 
the details of the optimization procedure, and section 4.2 discusses the distribution of the 
objective values for randomly generated controls. Section 14.31 computationally demon- 
strates the absence of traps in the control landscape for a five-level quantum system. Sec- 
tion U3] shows the dependence of optimization efficiency on the dimension of the quantum 
system, and Sec. 14 .51 examines the dependence of the computational complexity on the de- 
generacy structure of the eigenvalues of the initial state and target observable. Section I4~6l 
compares the computational efficiency of coherent and incoherent control. Optimization 
over constrained landscapes is investigated in Sec. Concluding remarks are given in 



2 Formulation of control for arbitrary TV-level sys- 
tems 

In this section the evolution of controlled iV-level open quantum systems is modelled 
by Kraus maps. As background, first the common formulation of the objective function 
in terms of Kraus operators is provided. Then the control problem is reformulated as 
optimization over a suitable Stiefel manifold; this representation is used in the subsequent 
numerical analysis. 

2.1 Kraus maps 

Let M.n be the vector space oi N x N complex matrices, with identity matrix 1^. The 
density matrix p G M.n of an iV-level quantum system is a positive semidefinite (and 
therefore Hermitian) matrix with unit trace, Trp = 1. A linear map $ : M.n — > M.n is 
positive if $(M) > for any M G M.n such that M > 0. The most general evolution 
transformations of density matrices are given by linear Kraus maps $ : J^-n — * -M-n-i 
which are defined by the following two properties: 

• Complete positivity: For any integer n, the map $ ® I n acting on M.^ <E> M. n is 
positive, where ® denotes the Kronecker product. 

• Trace preserving: For any MeMjv, Tr$(M) = TrM. 

Any Kraus map $ can be written in the Kraus operator-sum representation (OSR) 



Seed 



form 




(1) 



i=i 
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and the trace preservation condition implies for the Kraus operators Ki e M.^ that the 
relation is satisfied 

l 

Y,K\K t = l N . (2) 

i=l 

There exist many equivalent operator-sum representations of the same Kraus map. In 
particular, as shown in for any OSR with I > N 2 Kraus operators there exists an 
equivalent OSR with no more that N 2 Kraus operators. Thus, we only need to consider 
the OSR with I = N 2 Kraus operators (some of the Kraus operators can be zero matrices). 
Even for I = N 2 the decomposition ([1]) is not unique. Indeed, let lA{n) be the set of n x n 
unitary matrices, and let U G U(N 2 ) be a unitary matrix with matrix elements tiy. Define 
a new set of Kraus operators by the relation 

N 2 

K i = Y< u ^ 1<J<^ 2 - (3) 
i=i 

Then YliLi ^l^i = ^-N an d = £i!ti K-ipK-l = £$!ii K-ipK-l f° r an Y P- Therefore the 
two sets of Kraus operators {Ki} and {Ki} provide two equivalent representations of the 
same Kraus map. 



2.2 The objective function: formulation in terms of Kraus op- 
erators 

The optimization goal in quantum control is to maximize the objective function J = 
(0) = Tr[$(p)0] = Tr[p tf O], where is some target Hermitian operator, (•) denotes the 
expectation value at the final time t{, and pt { is the state of the system at the final time, 
evolved under controls from some initial state p = pt . The Kraus operators {Ki} = 
{Ki(tf, to)} describe the generally non-unitary evolution $ of the initial density matrix 
p at time t into a density matrix p tf at time tf, such that p tf = = ^^KipKj. 

They contain the information about the system-environment interaction, all control field 
interactions, and the state of the environment which also can be used as a control. Hence, 
J is a function of the Kraus operators 

N 2 

J[K X , . . . , K N i] = Tr [J2 KipKje] , (4) 

i=i 

and the control goal can be formulated as a constrained optimization problem: given p 
and O, maximize J over all sets of operators {Ki} that satisfy the constraint (T5]). 

For the remainder of the paper, we will take p and O to be simultaneously diag- 
onal. Indeed, we can always choose a basis in which O is diagonal, and write p and 
{Ki} in this basis. Since p is Hermitian, there exists a unitary matrix Q such that 
p = QcrQ\ where a is a diagonal matrix. Then the objective function (jl]) takes the form 

J = Tr £2=i KiaKjO , where K { = Kfr. The new Kraus operators {Ki} also satisfy 

the constraint ^ and the objective function is equivalently represented as a function of 
Ki with simultaneously diagonal matrices a and O. 
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2.3 The objective function: formulation in terms of Stiefel man- 
ifolds 

The above formulation can be expressed more succinctly in terms of the Stiefel manifold 
[10] ■ Let Ai(n, k,W) be the set of n x k matrices with matrix elements in the field F of 
real or complex numbers (i.e., F = R or F = C). The Stiefel manifold is defined as 

V k {¥ n ) = {Se M(n, k, F) : &S = I k }. 

The manifold \4(F n ) is called a real (resp., complex) Stiefel manifold if F = R (resp., F = 
C). Given a Kraus map $ and a set of Kraus operators {Ki}, we form the corresponding 
iV 3 x iV Stiefel matrix S as follows: 

( Kl \ 

S = K ! . (5) 
\K N J 

The constraint (T5]) can be expressed as the equality S^S = In, which defines the complex 
Stiefel manifold S = Va^C^ 3 ). Furthermore, the objective function can be written as 
a function of the Stiefel matrix S 

J(S)=Tr[SpS*(I N 2®G)], (6) 

The control goal in this formulation is to maximize the objective function ([6]) over the 
Stiefel manifold S. Note that the objective function (jSJ) is by construction real valued for 
any initial density matrix p and for any Hermitian target operator G. 

We now address the non-uniqueness of the Kraus operator parametrization in terms 
of the Stiefel manifold. Let W = {U <8> Ijv : U G U(N 2 )}. It is straightforward to verify 
that yS G S and VW G W holds S = WS G S. If {Kj} and {Ki} are two sets of Kraus 
operators that determine two Stiefel matrices S and S through (jSJ), then they define the 
same Kraus map and are related by the equality (j3J) if and only if 3W G W such that 
S = WS. Thus, equivalent parametrizations of the same Kraus map correspond to Stiefel 
matrices related by S = WS with some W G W. This property implies the invariance 
of the objective function under W-transformations, J(S) = J(WS) for any W G W, and 
will be used in Sec. [5] for analyzis of the search effort for optimization of J with additional 
constraints on the available Kraus operators. 

The Stiefel manifold Vfc(F n ) can also be defined as the set of orthonormal fc-frames 
in F n jH]. In this way, the Stiefel manifold S can be specified as the set of ordered 
iV-tuples Xi, . . . ,Xn G C^ 3 such that (X^Xj) = 5ij, where 5ij is the Kronecker delta 
symbol and (•, •) denotes the inner product in C^ 3 . In the remainder of the manuscript, 
the notation (•, •) will be used for inner products in several appropriate different spaces 
(namely, standard inner products in C^ 3 and in C^ 2 , and real Hilbert-Schmidt inner 
product in S and in the tangent space T S S at S). Vector Xj in this representation 
contains elements of the ith column of the Stiefel matrix (jSj) in certain order and can be 
decomposed in the direct sum 

X l = Yi®Yi®---® Y N 
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Here each Y- G , where 1 < i,j < N, is a complex vector of length N 2 of the form 
Yj = {(Ki)ji, (K 2 )ji, . . . , (K N 2)ji}, i.e., components of the vector Y- are the ji-th matrix 
elements of all the iV 2 Kraus operators ifj. The orthogonality condition (Xi,Xj) = 5ij 
implies the relation 

N 

fc=i 

Here (•, •) denotes the inner product in C^ 2 and should be distinguished from the same 
notations used above to denote the inner product in C^ 3 . 

The objective function for diagonal matrices p = J2^ =1 Pi\i) (i\ and 9 = ^2j =1 9j\j)(J\ 
can be written as 

N 

J[{Yj}} = £ Wii V*. (8) 

It is clear that ^min < J(S) < 9 max , where 9 m[n and 9 max are the minimum and maximum 
eigenvalues of 6, respectively. Indeed, we have 

N N 

9 mi n J2 W Y iW 2 Pi ^ J ^ E W Y ifPr 

i,j=l i,j=l 

Now, by first summing over i and using (jTj), and then summing over j and using Trp = 1, 
we have the desired inequalities. 

Since the maximal value of the objective function J equals to # max , the set of optimal 
controls (i.e., the set of all Stiefel matrices which maximize the objective function) is the 
manifold A4 max = {S E S : J(S) = 6* max }. For the case 9 = |iV)(iV| of special interest, it 
follows from flSJ) that 

-Mmax = {Oi>X=i : H 1 ?!! 2 = <W for any % such that Pl + 0} . 



3 Gradient and Hessian of J 

The numerical analysis in section 0] uses a gradient algorithm for optimization of the 
objective function J(S). This algorithm requires solving the equation 

^ = gradJ(S). (9) 

Here grad J is the gradient of the objective function, which induces the corresponding 
gradient flow on the Stiefel manifold S via Eq. (Q. 



3.1 Gradient of J 

We now derive an explicit expression for the gradient. Denote the differential of J at 
S G S by dsJ : TsS — > M, where TsS is the tangent space at S. By the product rule for 
derivatives 

d s J(SS) = KTr [{5S)pS\l N 2 ® 0) + Sp(5Sy (I N 2 g> 9)] (10) 
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where real part 3? is taken since the objective ([6]) is a real function. Since 31TyA = 3?TrA^ 
for any matrix A, the second term in the right hand side of fllOp can be rewritten as 
KTr[(I W 2 ® Q)(5S)pS^] and we get 



3?Tr [(^)p^t( Ijv2 g, 0) + ( Ijy2 g, 6)(55)pS t ] 
KTr [pS^I^ ® 6)( < J5) + pS^Ijv* ® 
23Wr [pS^Ijva <g> 6) (55)] 
(2(Ijv2 ®0)Sp, 55), 



(11) 



where SS G T^iS, and (A, B) := dtTr[A^B] is the inner product on S and TsS. By the 
Riesz Representation Theorem, there exists X G TgiS such that ds J(5S) = (X, 5S) for 
all SS G Ts«S. The vector X is the gradient of J at S 1 , denoted by grad J(S). 

Since grad J(S) must lie in TsS, it is necessary to remove the component orthogonal 
to TgS from the vector 2(1^2 <8> 0)Sp appearing in the last line of Eq. (fTTj) . Differentiation 
of the identity S^S = Ijy gives S^(SS) = —(8S)^S, so S*(5S) is skew- Hermit ian. This can 
be rewritten as SS = SB + (I N 3 — SS^)D, where B G M.(N, N, C) is a skew-Hermitian 
matrix and -D G .M(iV 3 , iV, C) is an arbitrary matrix. (Note that S^(SS) = B, since 
S^S = I N ). Any A G .M(iV 3 , N, C) can be decomposed as follows: 



Let C = {S^A + A^S)/2 and B = (S^A - A*S)/2, so that A = SC + SB + (I^a - Stf^A 
Clearly C is Hermitian and B is skew-Hermitian, so (SC,SB — (I^a — S'S'^A) = 0. 
Therefore, SC is orthogonal to TsS, and hence (A — SC) G Ts«S. As a result, 'Ps'(A) = 
A - S(S^A + A Jf S)/2 is an orthogonal projector from M(N 3 , N, C) onto T S S, and 

grad J{S) = 2(I N 2 Q)Sp - S [S\l N * <g> 6)5p + ((1^ <g> 0)5p) f 5] 



3.2 Hessian of J 

In the analysis thus far, we have only considered grad J(S), which gives first-order infor- 
mation about J. The Hessian gives useful second-order information about the minima, 
maxima, and saddles of J (where grad J(S) =0). At such points, the eigenvectors of 
the Hessian with positive (resp. negative) eigenvalues correspond to directions in which 
J increases (resp. decreases). 

The Hessian of J at S G S acting on 5S G TsS is defined as the covariant derivative 
of grad J(S) in the direction SS [12"] : 



Covariant differentiation of a function on a vector space is equivalent to taking the ordinary 
differential. However, S is not a vector space. In the following, the strategy will be to 
take the covariant derivative of grad J(S) as a function on M. {N 3 , N, C), which is a vector 




(2I N 3 - SS*)(I N 2 <g> Q)Sp - SpS\l N 2 <g> 0)S. 



Hess J(S) : T S S -> T S S, 



Hess J(S)(6S) = V^grad J(S). 
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space, and then project this onto S. Since S inherits its inner product from M(N 3 , N, C), 
this strategy gives the covariant derivative of grad J(S) on S. 

We now calculate an expression for the eigenvalues and eigenvectors of the Hessian of 
J on the critical manifolds. By differentiating grad J(S) in the direction of SS, we obtain 

V«?grad J(S) = d s grad J(5S) = 2(I N 2 <g> 0)(5S)p 

-{SS) S\l m <g> Q)Sp - S{5Sy {I N 2 <g> Q)Sp 
-SS\I N 2 <g> e)(5S)p - (8S)pS^(I N 2 <g> Q)S 

-Sp{5sy {i N 2 <g> 0)5 - sp&iiN* <g> e)(55), 

where V denotes the Riemannian connection on Ai(N 3 , N, C). We now project this onto 
Ts<S. Letting A = Vss grad J(S) gives 

Hess J(S)(6S) = V 5 sgrad J(S) = P 5 (V 55 grad 7(5)) = A - ^(SU + A^). 

With some algebra, this expression can be reduced to 

Hess J(S)(5S) = 2{I N 2 ® Q)(SS)p - (SS)S^(I N 2 <g> Q)Sp - (5S)pS\l N 2 <g> G)S 

+^[SS\5S)S^{I N 2 <g> Q)Sp + SS^{5S)pS^I N 2 <g> Q)S 

-2SS\l N 2 <g> e)(5S> - 2Sp{6S)i {I N 2 <g> 6)5 
+Sp&(I N * <g> Q)S(6S)*S + SS^(I N 2 <g> 0)^(55)^]. 

Combining the first two terms in the square brackets gives 

ss^s) [^(1^2 ® e)5p + ^(ijva <g> 6)5] 

= SS^S) [2^(1^2 <g> 0)Sp - S^gmdJ(S)] 
= 2SS\5S)S\I N 2 <g> 6)5p, 

since gradJ(S') = at a critical point. Combining the last two terms in the square 
brackets gives 

[S P S\I N 2 <g> Q)S + SS\I N 2 <g> 0) 1 Sp]((5 1 S) t ,S) 

= [2(Ijv2 <g> Q)Sp - grad 7(5)] (( < JS , ) t S') 
= 2{I N 2 <g> e)Sp(8S)iS. 

As a result, we have 

HessJ(S)(5S) = 2(1^2 ®9)(iS)p- (55)^(1^ ®G)Sp - (6S)p& (I N 2®0)S 
-SS\I N 2 <g> Q)(5S)p + SS^SS)^^ <g> G)Sp 
-Sp(5Sy (I N 2 <g> Q)S + {I N 2 ® Q)Sp{SSyS. 
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4 Numerical assessment of optimization efficiency for 
landscapes without constraints 

This section presents numerical simulations, including (a) an empirical demonstration of 
the absence of suboptimal traps in the control landscape, (b) an analysis of the dependence 
of optimization efficiency on the dimension of the system N, target operator B, and initial 
state p, and (c) a comparison between coherent and incoherent control. 

4.1 The optimization procedure 

We now describe the procedure for the numerical analysis of the controlled excursions 
over the landscapes without constraints on the controls. First, an adapted version of the 
algorithm in [13] is used to randomly generate an initial Stiefel matrix Sq with a uniform 
distribution on the Stiefel manifold S. After the initial Stiefel matrix is generated, the 
Runge-Kutta method built into MATLAB is used to solve Eq. (Q with the initial condition 
S(0) = Sq. The method relies on using a variable step size. The tolerances in the 
differential equation solver are set so that \\S^S — I/v|| < 2 x 1CT 4 at any given point in 
the trajectory. Integration is terminated when J(S) > (9 max — 0.01). 

The efficiency of the optimization procedure is measured by the two parameters: (1) 
the number r of a-steps taken by the differential equation solver in MATLAB to reach 
the objective value J > (0 max — 0.01) and (2) the path length A taken to get there. A 
higher number r of cr-steps corresponds to a more difficult optimization problem. Given 
the number r of cx-steps, the path length A is defined as 

T-l 

A = ^||5(z + l)-5(z)||, (12) 

i=0 

where || jS'H = a/ (S, S) is the norm on S. Similarly, a large value of A corresponds to a 
convoluted trajectory through S and indicates an inefficient optimization. 

To ensure statistical uniformity, for some simulations an average was performed over 
the initial state with a uniform distribution. Uniform sampling on the space of diagonal 
density matrices is implemented as follows. Let be the standard simplex, i.e., the 
set of all vectors z = (zi, . . . , z n ) G W 1 such that Zi > and Y17=i Zi = !• Xi = 
— log(dj) where is uniformly distributed on [0,1], so are exponentially distributed 
with parameter 1. Now let 

Vi = — X — , I <i <n. 

Xi + . . . + x n 

Then the random vector y — (yi, . . . , y n ) is uniformly distributed on the simplex [S] 
and the diagonal density matrix p with matrix elements pa = yi is uniformly distributed. 
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4.2 The statistical distribution of the objective for randomly 
generated controls 

In practical optimization of the objective function, either in the laboratory or through 
simulations with a numerical algorithm, the initial control is usually randomly generated. 
As the Stiefel matrices serve as the controls, we first analyze the distribution of the 
objective value for randomly generated initial Stiefel matrices. Fig. [1] shows the mean 
value Jo of the objective function J (So,p) = Tr[SopSQ(IlAr2 <g> 0)] for = |jV)(iV| as 
a function of the system dimension N for a uniform distribution of the initial Stiefel 
matrix So and uniform distribution of the initial diagonal density matrices p. For this 
case the mean value J equals to 1/N. To understand this result, let . . . , \N) be 
an orthonormal basis in the Hilbert space of the system such that \N) = |^) is the 
target state. The uniform generation of the Stiefel matrix So and initial density matrix 
p does not have a preferred state and thus preserves the symmetry between the states 
|1), . . . , \N). Therefore, in the final density matrix p' obtained by applying to p the Kraus 
map associated to Stiefel matrix So, the averaged (over uniform distributions of p and 

JV 

So) population pi of each of these states will be the same for all i. Since ^ Pi = 1 

i=i 

and pi = P2 = ... = Pn, we have Pi = 1/N for each i. Hence, the mean value of 
J = Tr[S pSl(l N 2 <g> 0)] = p N will be 1/N for = |^)(^| being a projector onto the 
target state |\&). 

Fig. [1] shows the decrease in the expected initial value of the objective function along 
with a decrease in the standard deviation with increasing system dimension N. Fig. [2] 
presents the detailed form of the distributions for the cases N = 2 and iV = 10, respec- 
tively shown in 2a and 2b, with a uniform distribution of So on the Stiefel manifold and a 
uniform distribution of p on the set of diagonal matrices. In this figure, the distributions 
of the values of the objective function J are produced using 10 4 randomly selected pairs of 
So and p. The results agree with the natural expectation that the efficiency of a randomly 
choosen control decreases with increasing complexity of the system. The figure also shows 
that as ./V rises the distribution of the objective values becomes more concentrated around 
the mean value. An open issue is to obtain an analytical expression for the distribution 
of the initial objective value Jq(Sq,p). 

4.3 Absence of suboptimal traps 

Let A be a topological space and / : A — > R. The function / is said to have a local 
maximum at xq G A if there exists an open neighborhood of xq, U Xq C A, such that 
Vx G U XQ , f(x) < f(xo) and yet there exists some x\ G A such that /(xi) > f(xo)- If 
A represents the space of all controls and / : A — > R is the objective function to be 
maximized on A, then a local maximum of / is called a suboptimal (or false) trap in the 
control landscape produced by /. 

The control landscape for the objective function J : S — > R defined by (jSJ) is known to 
have no traps [32] . Figure [3] numerically demonstrates this general fact for a particular five- 
level quantum system. In the figure, three different initial density matrices are considered: 
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a pure state p = |1)(1|, a randomly generated mixed state, and a completely mixed state 
p = I5/5. The control goal is to transform each of these states into the final state 
/9f = |5)(5|, which maximizes the expectation of the target operator = |5)(5|. As shown 
in the figure, in each case the gradient algorithm is able to find the control corresponding 
to the maximum value J = 1 of the objective function. The algorithm was not impeded by 
suboptimal local maxima, for their presence would have caused the algorithm to terminate 
at J < 1. Many other cases showed the same trap free behavior (not shown here). 

4.4 Dependence of the search effort on the dimension of the 
system 

We now analyze how the dimension N of the controlled system affects the optimization 
search effort. The goal is to numerically analyze the statistical dependence upon N of the 
number of steps r to reach convergence and the path length A. In this section, the target 
operator 9 = |iV)(iV| is the projector onto the state \N). To obtain reasonable statistics, 
for each N we average over 50 simulations of the optimization procedure with randomly 
(uniformly across S) generated S(0) and randomly generated initial density matrices p. 
We also analyze how the number of zero eigenvalues of p (henceforth denoted as do) affects 
the scaling of optimization efficiency with N. 

In the case of mixed p, we change p at the start of each individual simulation. Figures H] 
and plot r and A in two different ways in order to illustrate the issues driving the scaling 
efficiency. For each of the six curves in Fig. m the number of nonzero eigenvalues N — do of 
the initial state p remains fixed. Each curve labelled by iV — do corresponds, for example, 
to the control of a sequence of quantum systems prepared initially in a state at a relatively 
low temperature, with no population in do high eigenstates of the density matrix. Both 
r and A do not show any significant dependence upon N. It is clear that for fixed N, the 
search efficiency is greater for larger values of do- However, increasing do for fixed N — d 
does not change the slope of the curves in Fig. HI showing that the complexity of the 
search remains relatively insensitive to N. The most efficient control problem considered 
in the figure is the transformation of a pure initial state p = (j 7^ N) into a pure 

final state p f = \N)(N\. 

In contrast, for the simulations in Fig. \5\ do is held fixed for all N. This corresponds, 
for example, to the control of a sequence of quantum systems with the initial state at 
ever higher temperature as N rises, producing a large number of populated energy states. 
Both r and A increase quite sharply as N increases, as shown in Fig. [51 It is clear that 
for fixed N, the efficiency of optimization increases as do increases. However, as with 
Fig. HI increasing d does not change the slope of the curves in Fig. \5\ showing that the 
efficiency remains sensitive to N. The most inefficient search corresponds to the control 
goal of transforming a maximum entropy initial state with do = to a pure final state 
with do = N — 1, which agrees with simple intuition. 

The conclusion from Figures H] and [5] is that when G is a projector, the search efficiency 
decreases with increasing numbers of nonzero eigenvalues of p. The overall dimension of 
the quantum system N has little effect upon the search efficiency, provided that the 
number of nonzero eigenvalues of p remains fixed. The large standard deviations in both 
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figures are most likely caused by fluctuations of the initial Stiefel matrix S and of the 
parameters of the initial density matrix p not included in the number of zero eigenvalues 
d . 

The results in Fig.H]have practical relevance. In the laboratory a sequence of quantum 
systems with increasing N and a roughly fixed small number of populated energy levels can 
be arranged. Under these conditions, the results shown in Fig. Vindicate that search effort 
in the laboratory should not be very sensitive to the dimension of the quantum system 
under control. This behavior is generally consistent with the broad fingings that system 
and environmental complexity appear to have little effect on the number of iterations to 
reach successful control in the laboratory. 

4.5 Dependence of the search effort on the degeneracy structure 
of p and G 

In this section, we analyze the dependence of the optimization search effort on the degen- 
eracy structure of p and 0. Recall that M. max = {S G S : J(S) = # max }, where 9 max is 
the maximal eigenvalue of G. As shown in [32], the dimension of -M max is 

dim(.M max ) = 2(d + ei )iV 3 - (2d e 1 + 1)N 2 , (13) 

where do and e\ are the degeneracies of the zero eigenvalue of p and maximal eigenvalue 
#max of G, respectively. The dimension of the maximum manifold as a function of do 
and e\ is plotted on Fig. [6] (a). If d is close to N, then the initial state p is close to a 
pure state, and for e\ close to N, the target operator G is close to a constant multiple 
of the identity operator. Equation (TT3"j) and Fig. (a) show that large values of d and 
ei correspond to higher-dimensional maximum submanifolds (note that the dimension of 
the maximum manifold on vertical axis of Fig. [6] (a) increases in the downward direction). 

Figures [6] (b) and (c) show the dependence of efficiency of optimization upon do and 
e\. As do and e\ approach N, the efficiency of optimization increases rapidly. Comparison 
with Fig. [HI (a) shows a strong positive correlation between the dimension of the maximum 
manifold and the efficiency of optimization. This result is expected, since an increase in 
dim(A^max) corresponds to a larger target submanifold of optimal solutions. The presence 
of the positive correlation is illustrated in a more explicit way in Fig. [71 where the two 
parameters r and A characterizing the efficiency of optimization are plotted versus the 
dimension of the maximum submanifold. The dimension of the maximum manifold is 
determined by the pair (do, ex) and different pairs can produce the same dimension of the 
maximum manifold. Each point in Fig. [7J corresponds to a pair (do, ei). The figure shows 
the general trend that an increase in the dimension of the maximum manifold decreases 
the required optimization search effort; however the correlation is not perfect and different 
pairs (do, ei) and (d' , e[) with the same or almost the same dimensions of their respective 
maximum manifolds can have different values of the parameters r and A. 
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4.6 Comparison of coherent and incoherent control 

We now compare the efficiencies of coherent and incoherent control. The coherent control 
mechanism is implemented as follows. Let S(0) be defined by Ki = U/N for some 
U E U(N). It is shown in |X] that unitary Kraus maps form an invariant submanifold 
of S with respect to gradJ(S). That is, if S u = {S E S : SU E U{N) : K { = U/N}, 
then the solution to dS/da = gradJ(S) with the initial condition S(0) = S E S u will lie 
entirely in Su- Hence, solving the differential equation allows us to simulate density matrix 
evolution by coherent unitary control. Indeed, then p t = J2iLi Ki(f)pK\(f) = U{t)pU\t). 

In all the simulations here, = |iV)(iV|. For unitary control, the maximal value of 
the objective function J[U] = Tr[UpWQ] is the maximal eigenvalue p max of the initial 
state p. Thus, to ensure a fair comparison between the coherent and incoherent control, 
the target observable value is set to J = p max for both incoherent and coherent control, 
and the algorithm stops as soon as the value J = p max — 0.01 is attained. This stopping 
criteria is the reason for the difference between the curve corresponding to incoherent 
control in Fig. (a) and the curve in Fig. [HJ in the simulations displayed in the prior 
figure, the target observable value was J = 1 rather than J = p max . 

Fig. [8] shows that with the ability to generate arbitrary Kraus maps, incoherent con- 
trol can be a far more efficient process than coherent control for both pure and mixed 
p, especially for large values of N. The greater freedom allowed by incoherent control 
decreases the complexity of the problem and allows for a more efficient search. 

5 Control under linear constraints on the Kraus op- 
erators 

This section considers control under additional constraints on the available Kraus maps, 
which produce constraints on the Stiefel manifold. The target operator is assumed to 
have the form 9 = \N)(N\. 

Let h : S — > M. q be a set of q real-valued constraints. Recall from Sec. 12.31 that 
the objective function J is invariant under W-transformations. Since W-transformations 
correspond to different parametrizations of the same physical evolution Kraus map, any 
reasonable constraint should be W-invariant, and thus we impose the requirement that 
h(S) = h(WS) for any S E S and any W G W. 

We restrict the attention to affine constraints, which are of the form h(S) = g(S) —7, 
where g is linear over R and 7 e I 9 is a constant. Specifically, for a given set of matrices 
{Bi, . . . , B n } we consider W-invariant affine constraints of the form Tr(BjKj) = for each 

I = 1, . . . ,n, j = 1, . . . , N 2 , and B t E M N . Since Tr (BjKj) = Tr u^K^ = 

by linearity of the trace operation, this constraint is W-invariant and the set of Kraus 
matrices satisfying this constraint forms a W-invariant subset of the Stiefel manifold. 
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5.1 Numerical procedure 

The constraints Tr(BjKj) = with l<i<n, l<j< N 2 can be rewritten as a set of 
2nN 2 constraints h k : S -»• R defined as follows. Let Gi E M(N 3 , N, C) for I = n(j-l)+i 
be the matrix with Bi occupying rows (j — 1)N + 1 through jN and with other matrix 
elements set to zero. Then define G k = Gk for k = 1, . . . , nN 2 , Gk = iGk- n N 2 for 
k = nN 2 + 1, . . . , 2nN 2 and set hk(S) = (Gk, 5). The control goal is to maximize J over 
/' '«'). 

First, we need to find a matrix S E Sh which represents an initial control satisfying 
the constraint. To do this, define 

2nN 2 

f(S)=Y,(Gk,S) 2 . (14) 
fc=i 

We see that 

dsf(SS) =YZT 2 (Gk,6S)(G k ,S) = (j:TJi 2 2(G k ,S)G k ,8s). (15) 

Hence, grad/(5) = Y^k=i S)G k . Now generate an arbitrary So E S and solve the 

equation 

^ = -P ff (grad /(£)), 

with the initial condition S(0) = So, where Vs is the orthogonal projector from A4(N 3 , N, C) 
onto T$S (see Section \3.1\\ . Then, if the landscape of / on S is trap- free, the algorithm 
will always find a global minimum S of /, which will satisfy the constraint h(S) = 0. It 
is unknown whether this constrained landscape is trap free. 

After producing the initial Stiefel matrix S, we maximize the objective function J on 
Sh by solving the differential equation 

^ = 7V(gradJ(S)), 
dcr 

with the initial condition 5(0) = 5. Here grad J(S) is the gradient of J on S and Vh,s is 
a projector from T S S onto T s Sh- The expicit expression for Vh,s is derived in iBl 



5.2 Numerical results: general linear constraint 

It is difficult to derive a general analytical expression for the maximum value of J on 
the constrained manifold /i _1 (0) due to the complicated nature of the constraints. For 
this reason, we cannot determine that the gradient algorithm is stuck at a false trap 5 
(where grad J(S) = 0) by simply calculating J(S). Therefore, for a fixed constraint h, we 
performed the optimization procedure ten times using a different initial condition 5(0) 
each time and compared the resultant ten maximal values of the objective function. Let 
5* be the optimal control on the i th run (where 1 < i < 10), with corresponding maximum 
value J* = J(S*). If Jl < J* for some k and /, then SI is a false trap. Note that Jl = J* 
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for all k and I does not guarantee that the landscape is trap-free; the only conclusion is 
that the ten runs of the algorithm have not found a false trap. 

We performed simulations for N = 2,3, 4. For each N, five different initial states p were 
generated, and for each p five different collections of matrices {-Bi, . . . , B n } corresponding 
to five constraints were produced. We consider n < N 2 — N — 1, where iV 2 — N — 1 
represents the maximum number of constraints of special form corresponding to fixing 
to zero individual matrix elements of the Kraus operators. As a result of the numerical 
optimization, each of the ten runs performed with initial controls 5j(0) produced the 
same maximal value J* = J*, and therefore we did not find a false trap. Although this 
result does not prove the absence of false traps for linear constraints, it indicates that it 
is surprisingly difficult to find such traps, if they exist. 

5.3 Numerical results: fixing to zero individual matrix elements 

We now consider a special case of the W-invariant linear constraints such that h : S — > 
M? N is the constraint (K{)ij = for all I = 1,...,N 2 and for some pair The 
constraint corresponds to setting the ji th element in each of the A^ 2 Kraus operators to 
zero; we consider the real and imaginary parts separately, hence there are 2A^ 2 constraints. 
Since K n = J2 m=1 u nm K m defines the W-transformation, (K n )ij = for all n as well. 
Hence, h(WS) = h(S), and the constraint is W-invariant. More generally, we consider 
W-invariant constraints of the form 



(i^ M , = (O = 0, Z = l,...,iV 2 , V 3q eh, \/i q eI 2 , q = l,...,n, (16) 



where I\ and I2 two subsets of the set {1, 2, ... , N} each with n elements. 

For such a constraint, equation ([8]) can be used to determine analytically the optimal 
value of the objective function J on the constrained set /i _1 (0): 



For each A^ = 2, 3, 4, 5, we fix to zero n matrix elements of every Kraus operator, with n 
between A^ and A^ 2 — A^ — 1 (the maximum possible number of matrix elements which can 
simultaneously be fixed to zero). For a given n, the optimization procedure was performed 
25 times, and a different collection of matrix elements was fixed to zero during each run 
(i.e., different sets I\ and I2 were choosen). The gradient algorithm was able to reach 
the maximal value J max each time, showing that there do not appear to be false traps 
in this landscape. If suboptimal maxima were encountered, the algorithm would have 
gotten stuck at J < J max , and global optimization could not have been performed. Thus 
the optimization procedure did not discover any false traps for 25 randomly generated 
constraints. Again, this could not be taken as conclusive proof of the absence of false 
traps. Evidently, more complex or demanding constraints are called for to find traps. 



J, 



max 
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6 Conclusion 



This paper analyzes the efficiency of optimization over control landscapes for open quan- 
tum systems governed by Kraus map evolution. Several conclusions stem from the find- 
ings. When is a rank-one projector, which corresponds to the control goal of trans- 
forming an initial state p into a pure state, the search efficiency primarily depends on 
the number of nonzero eigenvalues of the initial state. The efficiency is relatively insensi- 
tive to the dimension of the quantum system N, provided that the number of populated 
energy states in the initial density matrix remains constant. As the number of nonzero 
eigenvalues of p rises with N, the search for an optimal control becomes less efficient. This 
result agrees with the expectation that transforming a high-entropy initial state into a 
low-entropy final state is a more difficult control problem than controlled transformations 
between states with similar entropy. 

The analysis also reveals that for fixed N, the search efficiency positively correlates 
with the number of zero eigenvalues of p. This result can be extended to a more general 
principle: when the dimension of the quantum system is fixed, the dimension of the 
maximum submanifold (the set of Kraus operators that correspond to optimal control) 
positively correlates with the efficiency of the optimization procedure. This statement 
agrees with the common intuition that a "larger" target results in an easier and more 
efficient search. The scaling behavior with iV found in this work is also consistent with 
that identified with unitary evolution, both dynamically and kinematically [35, 36J. 

We then showed that incoherent control modelled by Kraus map evolution, under 
the assumption that any Kraus map can be generated, is more efficient than coherent 
control modelled by unitary evolution. The larger number of control variables available 
in incoherent control actually decreases the complexity of the search effort. While the 
influence of the environment makes the total system ostensibly more complicated, the 
results show that the ability to control the environment can decrease the search effort. 

We also analyzed control landscapes with linear constraints on the Kraus maps. Even 
with the maximum possible number of linear constraints, false traps were not found. 
While this result does not prove the absence of false traps, it is nonetheless surprising. In 
the future work, we would like to investigate the control landscapes for constrained Kraus 
maps in more detail both numerically and theoretically. 

The kinematic analysis needs to be extended by a more detailed investigation of the 
role of the critical structure of the control landscapes on the search effort. In particular, 
the possible influence of saddle manifolds on the required search effort should be ana- 
lyzed. This analysis may reveal more subtle structural details about the quantum control 
landscapes. Also non-topological properties of quantum control landscapes may affect 
the optimization efficiency. In general, it is necessary to find all essential characteristics 
of the initial state p and the target operator that affect the efficiency of the search. 
Another important problem is to study the dynamics of controlled open quantum sys- 
tems with regard to topological and non-topological characteristics of the corresponding 
dynamical control landscapes. Various specific model systems can be used to study the 
dependence of search efficiency upon the parameters characterizing the system and en- 
vironment. The presence or absence of false traps in the dynamical control landscapes 
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should be investigated, including situations with constraints on the dynamical controls. 
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A Appendix A. Invariance of the submanifold Su for 

dS/da = grad J(S(a] 



Here we show that the submanifold Su := {Su G S \ 3U G U{N) such that Ki = j^U for i = 
1, . . . , A 2 } (i.e., all of the Kraus matrices determining a point Su G Su are equal to the 
same constant multiple of some unitary matrix) is invariant for the differential equation 
dS/da = giadJ(S(a)). 

Let A be a manifold with tangent bundle TX. Consider the differential equation 



dx/do = f(x(a)) 



(17) 



where / : X — > TX is a smooth function, and x : [0, 1] — > A is a path through A 
parametrized by the real variable a. A manifold Y C A is called an invariant submanifold 
for the differential equation (fTTj) if x(0) G Y implies that x(a) G Y for all a G [0, 1]. A 
compact manifold Y C A is an invariant submanifold for (1171) if and only if f(x) G T X Y 
for each x G F [15] . 

It was shown in Sec. 13.11 that 5S G T$S if and only if S^dS) is skew-Hermitian. 
Therefore, writing SS as a stack of A 2 A x A matrices 5 Si, . . . , 6S N 2, we see that for any 
Su G «Sw, SS G Ts v Su if the matrix [/t X^fcLi <^fc is skew-Hermitian. 

Theorem 1 Let E be the N 2 x 1 matrix (column vector) with all elements equal to one 
(i.e., E(z) = 1 /or oW ij. Then for any S = ^(E ® C7) G <S W; £/ie maim Z := ^[(E 1 " <g> 
Ijv) grad 7(5)] zs skew-Hermitian. 



Proof. Recall that grad J(S) = (21^3 

2 



S&Xlw <8> 6)Sp - SpS* {I N 2 <8> 6)5. Then 



[/ t (E t ® I 



AO 



1*3 - -^(E ® ?7)(E t <8> t/t) (I^ 2 ® 0)(E ® 
LA A d 



1 



-C/t — (E <g> C/)p(E t ® ^(I^a <g> 0)(E <g> C/) 
JM 6 



e)(E® c/> 



-Z7t-^(Et®£7t)(l^ 



e)(E® c/) 



|:C/ t (E t ® Ijv)Ijv3 - -^(E* ® f/ f )] (Ijva ® 0)(E (8) U)p 



0)(E® E7) 
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— (E f ® C/t)j ( E <g> 9C/)p - NpU ] QU 



= 2NU ] QUp - NU ] QUp - NpU^QU 

= N[U ] QU,p] (18) 

which is skew-Hermitian for Hermitian matrices p and 9. As a result, grad J(Su) G T Su Su 
for Sjj G 6^, so Su is an invariant submanifold for the differential equation dS/da = 
grad J(S(a)). 



B Appendix B. Derivation of the projector Vh,s 

Let h : {N 3 , N, C) — > M. q define a constraint on the Stiefel matrices, which restricts 
the set of addmissible controls to Sh . = S fl /i _1 (0). The goal is to find a projector 
7\s : TsS — > TsSh, such that the gradient of J on 5^ will be 7\s(grad J(S)). 
We will use the following lemma. 

Lemma 1 Lei X and Y be Riemannian manifolds and F : X — > Y '. Suppose that 
d x F is surjective for all x G X. Let fre £/ie operator on T X X defined as P x = I — 
(d x F)* o (d x Fo (d^F)*) -1 od x F. Then (a) P x zs a projection (that is, P x 2 = P x ) and (b) 

Proof, (a). It is straightforward to see that P% = P x . 

pi = (/-(d I F)*o(d/o(d/)r i od I F) 2 

= / - 2(d x F)* o (d x F o (d^F)*)- 1 o d x F 

+ (d x F)* o (d x F o (d.F)*)- 1 o d.F o (d x F)* o (d x F o (d^F)*)- 1 ° & X F 
= I- 2(d x F)* o (d x F o (d.F)*)" 1 o d x F + (d x FY o (d.F o (d^F)*)" 1 o d x F 
= I-(d x Fyo(d x Fo(d x F)*)- 1 °d x F = P x (19) 

(b). It is clear that if d x F(z) = 0, then z G T x F^ 1 (F(x)). Note that any vector v G T X X 
can be written as v = z + d^F*^), where w is arbitrary and z G T x F~ 1 {F{x)). Indeed, 
let w = (d x F o (d^F)*)" 1 o d x F(v). Then d x F(z) = d x F(v) - d x F o (d x F)* o (d x F o 
(d^F)*)- 1 o d x F(v) = 0, and therefore z G T x F~ l {F{x)). 

Now we will show that the image of P x lies in T X F {F{x)). For any v G T X X, write 
v = * + d x F*(w), where z G T x F^(F(x)). Then F r (f ) = F r (^) + P x o (d 2 F)*(w) = 
z - (d a F)* o (d,F o (d^F)*)- 1 o d x F(z) + (d 2 F)» - (d.F)* o {d x F o (d.F)*)" 1 o d.F o 
(d x F)» = z - (d x F)* o (d x F o (d.F)*)- 1 o d x F(z) = z since 2 G T^F^F^)) by 
assumption. Hence, the image of P x lies in T x F~ 1 (F(x)). This proves the lemma. 

Recall now that Vs is the projector from A4(N 3 , N, C) to T S S. Then dgh = dsh\ TsS 
and dg/i* = Vs d^/i*. If dgh is full-rank, then according to lemma [1] we have a projector 
Vh,s '■ TsS —>■ TsSh'- 

V h , s (SS) = dS-ds^oidshodshT'odsh^S) 

= 5S-V s o d s h* o (d s h oV s o dgh*)- 1 o d s h(5S). 
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In what follows, we will restrict our attention to affine maps denned by a set g = 
(gi, . . . , g q ) of bounded linear functional g^ : J\A(N 3 , N, C) — > R. By the Riesz Represen- 
tation Theorem, there exist unique matrices Gi G A4(N 3 , N, C) such that (ji(A) = (Gj, A) 
for all A E M(N 3 ,N,C). For a constraint of the form h k (S) = ^(Yj) k = (resp. 
3(Yj')fc = 0) as considered in Sec. 15.31 these matrices have the form G k = \d)(l\ (resp. 
G k = i\d)(l\), where d = (k - l)N 2 + (j - 1)N. 

Since 7 is constant and g is linear, dsh = g and dsh* = g*. To determine a formula 
for g*, note that for any i/Gl' 

q q q 

(g*(y),SS) = (y,g(5S)) = |>&(<fS) = ^(GUS) = (j^U^SS 

8=1 i=l i=l 

Therefore, g*(y) = Y^i=iViGi- Putting these expressions together gives 

(q \ / (GuEliViMGi)) 

d s hod s h*(y) = gl^yiPsiGi)] = \ 

Vi=1 ) \(G q ,EUyiMGi)) 

(G 1 ,V s (G l )) 



*=i \(G q ,Vs(Gi))J 
where Z has matrix elements = (Gi,Vs{Gj)). We finally get 

V h , s {SS) = 6S-V s og*{Z- l g(6S)) 
q 

= 5S-^'Ps(G i )(Z- 1 g(5S)) t 

1=1 

q 

= ^S-J2'Ps(G l )(Z-%~g,(6S) 
q 

= 8S ~Y,'Ps{G i ){Z-%{G i ,8S). 
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Figure 1: The mean value of the objective function J (S , p) = Tr[SopSl(I N 2 <8>0)], where 
Sq has a uniform distribution on the Stiefel manifold, p has a uniform distribution over 
the set of diagonal density matrices, and O = |jV)(iV|. 500 samples were taken for every 
point N. The error bars show the standard deviation for each N. 
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Figure 2: The distribution of the values of the objective function Jo (So, p) = 
Tr[S pSl(I N 2 <g> 0)] for iV = 2 [subplot (a)] and N = 10 [subplot (b)]. The initial Stiefel 
matrix So is uniformly distributed on the Stiefel manifold, p is uniformly distributed on 
the set of diagonal density matrices, and the target operator has the form = |iV)(iV|. 
10 4 samples were taken to produce the statistics for each plot. 
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Figure 3: The figure shows the value of the objective function J at each step o in the 
trajectory for a pure initial state p = diag(l, 0, 0, 0, 0), a mixed initial state, and a com- 
pletely mixed initial state p = 1/5 • diag(l, 1, 1, 1, 1). All three cases correspond to iV = 5 
and 6 = |5)(5| = diag(0, 0, 0, 0, 1). Each trajectory reaches perfect control J = 1 at the 
top of the landscape. 
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Figure 4: Dependence of search efficiency on the dimension of the system N for fixed 
numbers N — d of nonzero eigenvalues of p. Fifty simulations were performed for each 
point, and the average values of the number r of cr-steps and the path length A are plotted 
in (a) and (b), respectively. On each subplot, the six lines, from bottom to top, correspond 
to the number of nonzero eigenvalues of the initial density matrix N — d = 1, 2, 3, 4, 5, 
6, respectively. The error bars indicate the typical standard deviation of the data for the 
cases N — d = 1 and N — d Q = 6. 
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Figure 5: Dependence of search efficiency on N for different numbers d of zero eigenvalues 
of p. Fifty simulations were performed for each point, and the average values of the number 
t of a-steps and the path length A are plotted in (a) and (b), respectively. On each subplot, 
the eight lines, from top to bottom, correspond to the number of zero eigenvalues of the 
initial density matrix d = 0, 1, 2, 3, 4, 5, 6, N — 1, respectively. The error bars indicate 
the typical standard deviation of the data for the cases d = and d = N — 1. 
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Figure 6: Dependence of search efficiency on the degeneracy do of the zero eigenvalue of p 
and the degeneracy e x of the maximal eigenvalue of for N = 10. Figure (a) shows the 
dimension of the maximum submanifold as a function of these two parameters. Figure 
(b) shows the median number r of a-steps, and Figure (c) shows the median path length 
A. Figures (b) and (c) show that the search efficiency increases as the dimension of the 
maximum submanifold increases. 
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Figure 7: Dependence of search efficiency on the dimension of the maximum submanifold 
for N = 10. Figures (a) and (b) show the dependence of the median number r of cr-steps 
and median path length A, respectively, on the dimension of the maximum submanifold. 
The dimension of the maximum submanifold is determined by the pair (do, ei); each point 
on the plot corresponds to a (d ,ei) pair. These figures show the increase in the search 
efficiency as the dimension of the maximum submanifold increases. 
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Figure 8: Comparison of the optimization efficiency for incoherent control by Kraus maps 
and coherent control by unitary transformations, for both pure and mixed initial states p. 
Fifty simulations were performed to generate each point. The mean number r of cr-steps 
with typical standard deviations indicated by error bars is plotted on the vertical axis. 
The solid lines correspond to control by Kraus maps, and the dashed lines correspond to 
control by unitary maps. The lines marked by circles correspond to pure state p, and the 
unmarked lines correspond to mixed p. Similar behavior is observed for path length (not 
shown) . 



29 



